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ABSTRACT 

We present a method for simulating the evolution of HII regions driven by point sources of ionizing 
radiation in magnetohydrodynamic media, implemented in the three-dimensional Athena MHD code. 
We compare simulations using our algorithm to analytic solutions and show that the method passes 
rigorous tests of accuracy and convergence. The tests reveal several conditions that an ionizing 
radiation-hydrodynamic code must satisfy to reproduce analytic solutions. As a demonstration of our 
new method, we present the first three-dimensional, global simulation of an HII region expanding into a 
magnetized gas. The simulation shows that magnetic fields suppress sweeping up of gas perpendicular 
to magnetic field lines, leading to small density contrasts and extremely weak shocks at the leading 
edge of the HII region's expanding shell. 

Subject headings: HII regions — ISM: kinematics and dynamics — MHD — methods: numerical — 
radiative transfer 



1. INTRODUCTION 

Observations show that giant molecular clouds 
(GMCs) in local group galaxies convert at most a few 

S iercent of their mass into stars per cloud crossing time 
Zuckerman fc Evani Il974l ). and that clou ds are typ- 
ically destroyed in a few crossing times (Blitz ct ajj 
1200 9). well before they have converted a significant 
fraction of their mass into stars. Observations also 
strongly su pport the idea that GMC s are gravitation- 
ally bound (iKrumholz fc McKed [20051: iBlitz et all [20051: 
iKrumholz et al.l 120061 ) . so the fact that they survive 
for more than a crossing time yet do not collapse en- 
tirely into stars strongly suggests that internal feed- 
back plays a dominant role in CMC evolution. HII re- 
gions driven by newly-formed massive stars are likely to 
be the dominant sources of energy injection and mass 
loss i n clouds, and are therefore critical to GMC evolu- 
tion (iMcKcc fc Winiam ^'l997:'Wim ams fc McKe3 [l997l: 
lMatzner.,2002,) . IKrumholz et al. C200^ show using semi- 
analytic models that feedback from HII regions can quan- 
tiatively reproduce the observed lifetime, star formation 
rate, and star formation efficiency of GMCs. However, 
these calculations rely on simple analytic solutions for the 
evolution of spherically-symmetric HII regions in non- 
magnetic gas. Understanding the detailed evolution of 
GMCs will require a considerably more sophisticated nu- 
merical approach, and for this reason three-dimensional 
simulation of the evolution of molecular clouds under the 
influence of internal sources of ionizing radiation is a crit- 
ical problem in numerical astrophysics. 

In this paper we have the dual purpose of presenting 
a new algorithm for simulation of HII regions in molecu- 
lar clouds, and using this algorithm to explore potential 
computational problems that arise in general for simula- 
tions of this type. We also demonstrate our new algo- 
rithm in a simple application, the expansion of an HII 
region into a uniform neutral gas in which the magnetic 
pressure greatly exceeds the thermal pressure, as it does 
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in molecular clouds. In the past year several authors have 
considered the problem of simul ating HII regions, and 
presented both algorithms (e.g. [Arthur fc Hoard [20061 : 
[Mellema~et al. 2006) and results on the evolu tion of HII 
regions in turbulent hydr o dynamic media ([Dale et al.l 
[20051: iMellema et all [20051: [Mac Low et al.l [2006D . For 



a recent review see Hennev[ (|2006[ ). Several authors 
have also presented methods for simulation of ioniz- 
ing radiation- h ydrodyna mics (IRHD) i n a cosmological 
conte xt fe.g. [Abel fc Wandeltl 120021 : iWhalen fc Normal] 
[2006[) . which differs from the problem in the context of 
present-day molecular clouds primarily in the amount 
of cooling to which the gas is subjected and the condi- 
tions of the gas before it is ionized. Finally, researchers 
studying the evolution of ultracompact HII regions plan- 
etary nebulae have presented algorithms and results for 
ionizing radiative transfer with hydrodynamics and mag- 
netohydrodynamics in 2D and 3D under the simplifying 
assumption that the ionized gas has a perfectly sharp 
edge and is always in thermal and ionization equilibrium 
(e.g. [Garcia-Segura fc Franco|[199^[Garcia-Segura|[1997l : 
iGarcia-Segura fc L6ped[2000| ). 

Our algorithm improves on previous work in that it 
is the first to couple ionizing radiative transfer to mag- 
nctohydrodynamics (IRMHD) without imposing an as- 
sumption of thermal or ionization equilibrium. The in- 
clusion of magnetic fields is important because obser- 
vations indicate that the magnetic energy in molecular 
clouds is compa rable to the kinetic and gravitational po- 
tentia l energies (Crutcher 1999. 2005': [Heiles fc Crutcheil 
[2005() . Thus, dynamical expansion of ionized regions may 
be significantly altered by magnetic confinement, an ef- 
fect we wish to explore. Indeed, we show that even in 
the simple case of expansion of an HII region into a uni- 
form, magnetized medium, the magnetic field produces 
qualitatively new phenomena. The Alfven speeds of a 
few km s^^ typically found in molecular clouds reduce 
the strength of shocks associated with expanding ioniza- 
tion fronts at early times, and at later times turn off 
the shocks entirely, greatly reducing the collection of gas 
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and possibly thereby reducing the amount of triggered 
star formation. A non-equihbrium treatment of the ion- 
ization structure is important because we find that a 
non-equihbrium treatment of the thermal and ionization 
structure of the ionized gas leads to modest but signifi- 
cant effects on quantities such as the expansion rate of 
HII regions, effects that are lost in the assumption of 
perfect radiative equilibrium. 

Before exploring these new phenomena, however, we 
point out that there has yet to be a detailed study of 
potential computational problems and constraints that 
arise in 3D simulations of IRHD and IRMHD with the 
strong cooling and large temperature contrasts that are 
expected for modern day (as opposed to primordial) in- 
terstellar chemistry. We have therefore performed a de- 
tailed comparison of simulations using our method to 
analytic solutions for computationally challenging prob- 
lems, as a way of searching for potential difficulties that 
may arise in general IRHD methods. We discover several 
conditions that a simulation must satisfy to reproduce 
analytic results correctly. One must limit the amount by 
which the gas pressure is allowed to change between hy- 
drodynamic updates, for a ray-tracing method one must 
periodically rotate the orientation of the rays, and one 
must either resolve the ionization front or restrict the rate 
of cooling at the front to suppress excess cooling due to 
numerical mixing. We show that failure to meet these 
conditions produces quantitatively incorrect results. 

The remainder of this paper proceeds as follows. In 
§ [2] we describe the physical formulation of the prob- 
lem that we adopt, including our approximations and as- 
sumptions. In §[3] we present our simulation algorithm. 
In § 2] we compare our code to analytic solutions, and 
thereby demonstrate the existence of conditions that nu- 
merical methods must satisfy in order to reproduce ana- 
lytic solutions correctly. We use our method to simulate 
the evolution of HII regions in magnetized media in § [5l 
and finally we summarize and present conclusions in § [Bl 

2. PHYSICAL FORMULATION 

We wish to simulate the evolution of a magnetized 
molecular gas subjected to an ionizing radiation field due 
to point sources within it. The n sources of ionization are 
located at positions x„, and each has an ionizing lumi- 
nosity of s„, in units of photons per unit time above the 
Lyman limit. We treat the ionizing sources as monochro- 
matic, and assume that they emit all their ionizing pho- 
tons at frequencies near the ionization threshhold of hy- 
drogen. Since our focus is on the dynamics of the cloud 
rather than the detailed chemistry of the ionized gas, this 
approximation is quite reasonable. 

We describe the gas by a total mass density p, density 
of neutral species p„, velocity v, and total energy density 
E. The gas is threaded by a magnetic field B. The 
interface between an HII region and a molecular cloud is 
a photodissociation region (PDR) , within which there is a 
significant amount of atomic gas. In principle we should 
therefore separately track ionized, atomic, and molecular 
hydrogen, and possibly other species as well. However, 
PDRs do not have a signicant effect on the large-scale 
dynamics of the ionized gas unless the ionizing source is 
very weak or the neutral gas is much less dense than as 
typical for molecular clouds, as we show in § 12.41 For 
this reason, we neglect atomic gas, and assume that the 



density of molecular gas is p„ and the density of ionized 
gas is p- pn- 

2.1. Magnetohydrodynamics 

We compute the evolution of our gas using the equa- 
tions of multi-species ideal MHD including radiative 
heating and cooling and chemical evolution terms. In 
conservative form, these are 



| + V.(pv).0 
^(pv) + V • (pvv - BB) + VP* = 

dm 

— + V • (vB - Bv) = 

— +V.[{E + P*)v - B(B • v)] =g - £ 

^ + v•(p„v)=7^-x, 



(1) 

(2) 
(3) 
(4) 
(5) 



where P* = P -I- (B • B)/2 is the total pressure, P is the 
gas thermal pressure, TZ and X are the rates of recombi- 
nation and ionization, measured in mass per unit time 
per unit volume, and Q and C are the rates of radiative 
heating and cooling, in energy per unit time per unit vol- 
ume. The total energy density E not including atomic 
or molecular binding energies is 



E : 



B B 



(6) 



where e is the gas thermal energy density. We adopt 
an ideal gas equation of state, so the pressure is P = 
(7 — l)e with 7 = 5/3, corresponding to a monatomic 
gas. This choice of 7 is appropriate because ionized gas 
is monatomic, and hydrogen in molecular clouds, while 
diatomic, is generally too cool to access its rotational or 
vibrational degrees of freedom. It therefore acts as if it 
were monatomic for the purpose of dynamics. In addition 
to the evolution equations, the magnetic field is subject 
to the divergence-free constraint 



V • B = 0. 



(7) 



Physically, equations ([T]), ([2]), and (|4]) represent conser- 
vation of mass, momentum, and energy, equation ([3]) 
represents the approximation that the magnetic field is 
perfectly frozen in to the fluid, and equation ([5]) states 
that the mass of neutral gas is conserved by advection, 
and is altered only by ionizations and recombinations. 
Equation ([7]) expresses the non-existence of magnetic 
monopoles. Note that we have adopted a system of units 
in which the magnetic permeability fj, = 1; to convert to 
cgs units, we multiply our magnetic field strengths by a 
factor of v/Jtt. 



2.2. Recombination and Ionization 

To complete the evolution equations, we must spec- 
ify the rates of the radiative processes: recombination, 
ionization, heating, and cooling. Recombination is the 
simplest of these. We adopt the on-the-spot approxima- 
tion whereby any recombination to the ground state of 
hydrogen is assumed to yield an ionizing photon with a 
very short mean free path, which will in turn cause an- 
other ionization at essentially the same location as the 
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recombination (| Ost erbrockl 1 989f) . On the other hand, re- 
combinations to excited states of hydrogen produce pho- 
tons with long mean free paths, which are assumed to 
escape the HII region. Thus, the recombination rate is 
approximately 



7^ 



(8) 



where fj,H ~ 2.34 x 10~^^ g is the mean gas mass per 
hydrogen atom (for a gas of hydrogen and helium in 
the standard cosmic abundance), Ue and are the 

number densities of electrons and hydrogen nuclei, and 
a(^) is the case B recombination coefficient, ' ~ 
2.59 x lQ-^^iT/lO'^ ^"^ (|Osterbrocklfl989l : 

iRiikhorst et al.|[2005D at a gas temperature T. The H+ 
and e~ number densities are — {p ~ Pn)/ ^J^H and 

rie = +pac /{l4:fiH), where ac » 3 x 10~^ is t he car- 
bon abundance in the ISM (Sof ia fc Me"ve^l200lD . This 
assumes that all carbon is singly-ionized. Since the first 
ionization potential of carbon is smaller than that of hy- 
drogen, this is likely to be the case in regions where there 
is any ionized hydrogen present, which are the only re- 
gions for which we care about the recombination rate. 

In the on-the-spot approximation, the ionization rate 
is relatively easy to compute as well, since we need only 
solve the radiative transfer equation only along rays from 
ionizing sources. The ionizing flux at position x due to 
the nth ionizing source is 



F„ 



47r |x — x„ 



-r (x,Xt, 



where 



t(x,x„) 



{a-riH + cTj^n) d£ 



(9) 



(10) 



is the optical depth to ionizing photons along the path 
from x„ to X, nn = Pn/fJ-H is the number density of neu- 
tral hydrogen atoms, n = p/p-H is the number density of 
hydrogen nuclei, a = 6.3 x 10^^* cm^ is the cross section 
for absorption of a photon at the ionization threshhold 
by a neutral hydrogen atom, and 0-^ w 6.1 x 10^^^ cm^^ 
per H atom is the cross secti on for absorption of ion - 
izing photons by dust grains (jBertoldi fc Draind Il996l ) . 
For the purposes of the tests we report in this paper, for 
which we are comparing to analytic solutions derived for 
pure hydrogen, we set ad — 0- The photoionization rate 
is just the rate at which these photons are absorbed, so 
the total photoionization rate due to all sources at posi- 
tion X is 



X, 



ph 



-r (x,x„ 



47r |x - 



(11) 



We also compute the rate of ionizations due to collisions 
between neutral hydrogen atoms and electrons. 



2^coU ~ kcollTlePm 



(12) 



where the collisional io nization rate coefficient is 
(jTenorio-Tagle et al.lfT986[ ) 

fccoii ~ 5.84 X lO-ii^^e^-^-e eV)/(fcBT) ^^3 ^-i^ 

although for the gas densities and temperatures in the 
systems with which we are concerned here this term is 



generally negligible. The total ionization rate is the sum 
of these two terms. 



1 = Pn 



47r |x — Xr, 



(14) 



It should be noted that iRitzerveldl ()2005[ ) has recently 
pointed out that diffuse photons (i.e. those that have 
been absorbed and re-emitted) can be as important as di- 
rect photons in determining the structure of HII regions, 
so the on-the-spot approximation we adopt is potentially 
problematic. However, the importance of diffuse versus 
direct photons depends strongly on the degree of central 
concentration of the gas inside the HII region. Once an 
HII region becomes D type, which is the phase in which 
we are most interested from the standpoint of radiation 
hydrodynamic evolution, the gas inside the ion ized region 
is_nearly uniform in density. For this case, IRitzerveldl 
finds that direct photons dominate diffuse ones 
over the majority of the HII region volume, and the on- 
the-spot approximation is reasonable. Nonetheless, the 
on-the-spot approximation does make the shadows cast 
by dense gas sharper than they should be in reality, and 
we must keep this limitation in mind. 

2.3. Heating and Cooling 

We can approximately determine the heating and cool- 
ing rates by breaking them up into heating and cooling 
associated with ionization and recombination of hydro- 
gen atoms, and all other sources of heating and cooling. 
The heating rates due to ionization and recombination of 
hydrogen are quite simple in our monochromatic treat- 
ment of the ionizing radiation. In this approximation, 
each absorption of an ionizing photon deliv ers er ~ 2.4 
eV o f thermal energy to the gas (jWhalen fc NormanI 
l2006( l. so the photoionization heating rate is 



Sph = erCTUH 



E 



47r |x — x„ 



^-t(x,x„) 



(15) 



Each recombination of a hydrogen atom allows pho- 
tons to escape, leading to a loss of energy. A micro- 
physical calculation of th e resulting coolin g rate gives 
where (|Osterbrocklll989D 



Arec « (6.1 X 10-1" cm3 s-i) ksT 



-0.89 



(16) 



for temperatures T ^ 100 K, and recombination cooling 
falls to negligible values at lower temperatures. 

For heating and cooling that are not directly due 
to ionization of hydrogen, we use simple optically-thin 
heating and cooling curves. In molecular gas, we 
adopt the a pproxima te coo ling and heating functions of 
iKovama fc Inutsukal (|2002f) . and in p artially ionized ga s 
we compute the cooling rate foUowing lOsterbrockl (|l989l ). 
Our calculation includes the cooling by ion-electron col- 
lisions involving the first and second ionized states of O, 
N, and Ne, which are the dominant coolants in HII re- 
gions under solar metallicity conditions. We also include 
free-free cooling. Thus, the total heating and cooling 
rates are 



Q = erauH 



^4^|x-x„|' 



e"-(x.x„) ^ ,^^rKi (17) 



4 



Krumholz, Stone, & Gardiner 



and 

= AKi{T)njj+A,cc{T)nenH++J^ion-s{T)nenH+ (18) 

where the subscripts KI, r ec, and ion-ff indicat e heat- 
ing and cooUng following the lKovama fc Inutsukal curves, 
from recombinations, and from ion- neutral collisions and 
free-free emission in ionized gas. With these cooling 
curves, the equilibrium temperature in neutral gas at a 
density of nn — 100 cm~^ is approximately 11 K, and 
the equlibrium temperature in fully ionized gas is approx- 
imately 6400 K, with some weak density dependence. 

2.4. Atomic Gas and PDRs 

Here we justify our earlier statement that we can ne- 
glect PDRs in studying the large-scale dynamics of most 
HIT regions. We do this u sing the anal ytic PDR models of 
iBertoldi fc Draind (|1996[ ). Note that [Bertoldi fc Draind 
consider the case of convex fronts, in which a neutral 
cloud is bathed in ionizing radiation, whereas we are 
more concerned with concave fronts, for which the source 
is embedded within a molecular medium. To avoid com- 
plications arising from this, we will only consider the case 
of plane-parallel ionization fronts. However, the results 
should qualitatively generalize to our concave case. 

First consider a system a short time after an ionizing 
source turns on, before a shock front forms. In this case, 
there is only a dissociation front and an ionization front. 
The propogation speed of each front is then determined 
solely by the rate at which ionizing or dissociating pho- 
tons reach it. Let SLy and spuv be the stellar luminosi- 
ties in units of ionizing Lyman continuum (A < 912 A) 
and dissociating FUV (912 A < A < lllOA) photons per 
second emitted by the ionizing source, and TLy and rpuv 
be the optical depths from the source to the respective 
fronts. In this case an ionization front at a distance 
from the source propogates at a velocity 

w,: = , (19) 



where tih is the number density of hydrogen nuclei out- 
side the front. The dissociation front propogates at a 
rate 



Wd«0.15 



spuve 



47rr2n// 



(20) 



where is the front radius, and the factor of 0.15 arises 
because only a minority of FUV photons actually pro- 
duce dissociations. For stars of type 03 or earlier, i.e. 
those that drive significant HII regions, SLy/ sfuv ~ 1—2, 
so Vi > Vd unless the Lyman continuum optical depth is 
large enough compared to the FUV optical depth so that 
the fraction of Lyman continuum photons reaching the 
front is « 8 — 15% of the number of FUV photons reach- 
ing the front. This requires TLy — rpuv ^ 1-9 — 2.6. 

Since the dust extinction cross section for FUV pho- 
tons is larger th an for Lyman continuum photons by a 
factor of ~ 1.25 (|Bertoldi fc Drainill996( ). this can only 
happen if there is significant attenuation of Lyman con- 
tinuum photons due to recombining gas within the HII 
region. At early times before recombinations are signif- 
icant compared to photoionizations, this means Vi > Vd 
and the ionization and dissocation fronts will be coinci- 
dent. In this case, our neglect of the dissocation front is 
obviously not a problem. 



As the ionization front radius approaches the 
Stromgren radius. 



f 3sLy 



1/3 



\Awa(B)n%J ' ^^^^ 

recombinations within the ionized region produce a sig- 
nificant neutral population that raises TLy but not rpuv 
and slows expansion of the ionization front relative to the 
dissociation front. If this process continued indefinitely 
then the dissocation front would eventually race ahead of 
the ionization front. However, slowing of the ionization 
front will also eventually allow a shock to form, convert- 
ing it from R type to D type, and this will sweep gas up 
into a dense shell and modify the propagation speed of 
both the ionization and dissociation fronts. How much 
can the ionization front slow down before a shock forms? 
The condition for a shock to form is that the ionization 
front speed drops to of order the speed of a strong shock 
driven by the ionized gas, which is twice the ionized gas 
sound speed q. Thus, a shock forms when 



SLye" 



2& 



(22) 



If we now approximate that w once TLy 1 , it fol- 
lows that 



1 / SLyOf^^' HH 



I 2 ( 

: 4.1 In I S49n2Q!4 



(23) 



(24) 



where in the second step we have scaled to typical 
values for Galactic molecular clouds and HII regions, 

S49 = SLy/(1049 S-1), 712 = "ff/(100 Cm-^), ) ^ 

"^^Vajt'io'' K = a^-^V(2.59 x lO-^^ cm-^ g-i), and 
C6 — Ci/(10^ cm s~^). Thus, the dissociation front will 
not begin moving faster than the ionization front until 

TLy 1.9 — 2.6, but a shock will form once TLy ~ 4. 

Once a shock forms, wc must consider how the width 
of the PDR compares to the width of the shocked layer 
to determine if the dissociation front can ever pro- 
pogate ahead of the shock. In equilibrium, the column 
density of hydrogen ato ms through a PDR is roughly 
(|Bertoldi fc Draind[l996l ) 

^//«^?uv[l + (2-7 + 0o')*/']-\ (25) 

where ctfuv ~ 7.6 x 10~^^ cm^^ is the attenuation cross 
section of dusty gas to FUV photons and the factor 0o is 
a dimensionless number describing the balance between 
H2 formation and dissociation. Since the numerical fac- 
tor involving (^q is at most 0.21, this implies that there is 
a maximum hydrogen column iVpoR ~ 2.8 x lO'^^ through 
a PDR. Non-equilibrium PDRs are generally narrower 
than this. The column of hydrogen swept up in the dense 
shell around an ionized region shortly after it forms is 
roughly hhTs- If we compare this to the maximum col- 
umn density of a PDR, we find that the column of hy- 
drogen nHTs present when the shock first forms exceeds 
the maximum column of a PDR whenever 



SLy 



""3 uh 
= 2.3 X lO^^af 



(26) 
(27) 
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Thus, for Shy ^ 10^^ s~^, which is the case for any strong 
ionizing source, the dissociation front will essentially al- 
ways remain trapped between the ionization front and 
the shock front, except perhaps for a very brief interval 
when the front is changing from R type to D type and 
the Lyman continuum optical depth through the ionized 
region is about 2 — 4. As a result, the dissociation front 
will not have a significant effect on the dynamics of the 
ionized region, since the only gas it can affect will be 
gas that has already been shocked and is about to be 
ionized. The dissociation front will simply be passively 
carried along. 

This can break down if we consider weak ionizing 
sources, or if we have a medium substantially less dense 
than TiH = 100 cm~'^, the typical mean density in a 
GMC. If either of these is the case, then a dissocation 
front may run ahead of the shock front and pre-heat the 
gas, thereby affecting the dynamics. In this case our ne- 
glect of the atomic gas becomes more problematic. How- 
ever, even when pre-heating occurs, the dissociation front 
will only precede the shock front for a limited time, until 
the ionized region has swept up the critical column den- 
sity A^pDR- Thus, even for weak ionizing sources or low 
density regions, our approximation will only be invalid 
around the time when the front transitions from R type 
to D type. 

Our simple analytic calculation is in good agree- 
ment with the detailed one- dimensional simulations of 
iHosokawa fc Inutsukal ()2005[ ). who also find that the 
molecular hydrogen dissociation front around an expand- 
ing HII region becomes trapped between the shock and 
ionization fronts shortly after the front transitions from 
R type to D type. 

3. SIMULATION ALGORITHM 

3.1. Program Flow 

We solve the evolution equations on a Cartesian grid 
with uniform cell spacing Ax using an operator split ap- 
proach, in which we alternate conservative MHD updates 
with the source terms on the right hand sides set to zero 
with radiation updates in which we add the source terms. 
In designing this update algorithm, we bear in mind two 
factors. First, due to the comparitively simple nature 
of the radiation update, computing a single explicit up- 
date using the radiation source terms is computationally 
much cheaper than computing a single magnetohyrody- 
namic update. For this reason, we choose to sub-cycle 
the radiation update relative to the MHD update. This 
effectively makes our radiation update step implicit, since 
during any given simulation time step we iterate the tem- 
perature and chemical state of each cell over many ex- 
plicit update cycles. 

Second, the overall update time step is restricted by 
both radiation and MHD. The time step must obey the 
Courant condition, but we also wish to impose a heat- 
ing/cooling time step condition. If the energy in a cell 
changes too much between MHD updates, the corre- 
sponding large changes in pressure from one time step 
to the next may lead to inaccurate solutions. We there- 
fore wish to limit the amount by which the energy in a 
cell can change between MHD updates, which in turn 
imposes a time step constraint associated with the rate 
of heating or cooling. Note that, in contrast, the MHD 
update is not directly affected by the neutral density p„. 



so we need not impose a corresponding constraint on the 
overall simulation time step arising from changes in the 
chemical state of the gas. 

In order to satisfy these constraints, we perform a two- 
step update procedure, illustrated schematically in Fig- 
ure [T] At the end of time step n, in each cell we have a 
vector of quantities 

p" 

q" = B" (28) 

\Pn) 

at time t". (Here p is the gas momentum, and we are 
glossing over the fact that in the Athena algorithm the 
magnetic field is actually stored at cell faces rather than 
cell centers.) We wish to determine the properties q"+-'^ 
at time We first update the gas energy E and 

neutral density p„ using only the terms on the right hand 
sides of equations lU - ©. This update covers a time 
Ai" — min(At(7, Ate), where Ate is the Courant time 
step for the current configuration, and At^ is a time step 
imposed by a requirement that the internal energy of a 
cell not change by too much per time step. This gives 
us an intermediate vector q"'*. Since we do not want 
to have to perform an MHD update every time we do 
a radiation update, the radiation update procedure is 
iterative. Therefore E"-'* is not simply E" + At"- (G - £)"■ , 
and similarly for p^'* . We give a step- by-step description 
of the radiation update procedure in § 13. 2[ and defer 
formal definition of Ai^, and p^'* until then. 

Second, we perform an MHD update to gas quanti- 
ties starting from state q"'*, advancing the conservation 
equations ([I]) - ([1]) with the source terms on the right 
hand sides set to zero. This advance covers the same 
time At" as the radiation update. We perform this up- 
date using the standard Athe na explicit MHD update 
proce dure, which is described in lGardiner fc Stoni (|2005l 
l2006l ) , with the trivial addition of the passive scalar ad- 
vection equation ([5|) for the neutral density. The overall 
program flow is independent of the details of the MHD 
update algorithm. For the runs we report in this pa- 
per, we configure Athena to use its directionally-unsplit 
6-solve upwind constrained-transport method, the Roe 
Riemann solver with H correction, and second-order spa- 
tial accuracy. We do not describe these methods in 
greater detail here, and instead refer readers to the Gar- 
diner & Stone papers. At the end of this step, we have 
updated q" for all the terms in the evolution equations 
(HI) - ([5]), so the final state is q""'"^, and we may begin 
the next time step. 

It should be note d that our update cycle is quite similar 
to that of Whalen fc NormanI ()2006l ). 

3.2. Radiation Update Method 

Our radiation update algorithm combines ele- 
ments of the appr o aches lAbel fc Wandeltl ()2002D and 
IWhalen fc NormanI ()2006f l. together with some novel 
techniques that through the tests we describe in § |4] we 
found necessary to achieve high accuracy. Our procedure 
to advance the input state q" to the intermediate state 
q"'*is as follows. 

Wc do two prepatory steps before beginning a radia- 
tion update. First, we compute the Courant time step 
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Atp for the configuration q". Second , we construct a 
tree of rays following the prescription of lAbel fc WandelB 
(|2002i ) for each radiation source. We will not review 
the full formalism of this method here, but a summary 
of its relevant properties is that from each source one 
sends out a tree of rays in directions that correspond to 
the centers of pixels in the HealPix pixelization scheme 
(|G6rski et al.| [2005). At its coarsest level the tree con- 
tains 12 rays uniformly spaced over the unit sphere, but 
the scheme is hierarchical, such that any ray may be sub- 
divided into four child rays an arbitrary number of times. 
After the Zth such division, the unit sphere is discretized 
into 12(4') pixels, so that each pixel covers a solid angle 
J7 = (7r/3)4~' sr. As we follow rays outward, we sub- 
divide them to guarantee that at a distance r from the 
source fl < (l/2)Ax^/(47rr^), i.e. that the sohd angle 
that a given ray represents is always smaller than the 
solid angle subtended by a cell of the Cartesian grid by 
at least a factor of 2. We always divide each ray at least 
twice, so we never have fewer than 192 rays. 

Along each ray we construct a list of the cells through 
which that ray passes, and the length of the ray segment 
that lies within each cell. If we are running a calcu- 
lation in parallel, so that the computational domain is 
subdivided into a series of processor domains that are 
contiguous rectangular blocks, we first construct the list 



of processor domains through which a ray passes. On 
each processor we only store information about rays for 
which either the ray, its parent, or one of its children in- 
tersects that processor's domain. This enables us to pass 
information from one processor to the next as we walk 
outward along the ray tree without storing the full tree 
on each processor. 

If the source is moving, we repeat this procedure every 
time the source moves, generally every time we perform 
an MHD update. For a source that is fixed in space, we 
build a new tree every fifth MHD time step. Each time 
we build the tree we rotate it into a random orientation 
to minimize errors due to discretization in angle. We 
discuss this procedure, and why it is necessary, in more 
detail in § [43:31 

With these preparations done, we begin the iterative 
update procedure. We compute the ionizing flux passing 
through every cell by walking ou tward along the ra ys, fol- 
lowing the procedure outhned in lAbel fc Wandefl (2002) . 
Suppose a particular ray representing a solid angle Q 
starts at a distance tq from the source, and that the flux 
of ionizing photons at this point is Fq. The ray passes 
through a sequence of / cells, and we define £i as the 
length of the ray segment that intersects the ith cell, 
Ti — To + X] }=o the distance from the source to cell 
i, and Fi as the radiation flux reaching cell i. The op- 
tical depth through cell i is = {crnH + <Tdn)ii, where 
riH and n are the number densities of neutral hydrogen 
atoms and hydrogen nuclei in the cell. The flux reaching 
the next cell is just 



- F,e- 



'i+i 



(29) 



and conservation of ionizing photons therefore implies 
that the number of ionizations per unit time in cell i 
due to photons moving along this ray is Fi{l — e^'^')flrf . 
The corresponding rates of mass photoionization per unit 
volume and photoionization heating per unit volume are 

2 



Jph^F,il-e-^^)nr 

gpb^F,{i~e-^^)nr 



,2 er 



(30) 
(31) 



We use this procedure to trace the photons along each 
ray, starting from the 12 coarsest ones. Once we have 
reached the end of a ray, we pass the flux Fj that escapes 
the flnal cell on to its child rays and repeat this procedure 
for them. We continue until we either reach the edge of 
the computational domain or until such a small fraction 
of the photons remain (we use 10"'^) that it is no longer 
worthwhile to continue following them. The total mass 
photoionization rate Tp^™ and photoionization heating 
rate t/pj^™ iii ^ cell at time i"'™ is simply the sum of 
the rates contributed by all the rays that intersect it, 
as determined by equations (|30p and (PT|) . for the gas 
properties at that time q"'™. 

Once we have determined the rates of change of the 
neutral mass and energy due to photoionization, we de- 
termine the rates due to recombination, coUisional ioniza- 
tion, and optically thin heating and cooling using equa- 
tions (HI), (|T2l), dlel), dni), and These calculations 
are purely local, and we simply perform them using the 
current state of the gas q"'™. For reasons we discuss in 
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§ 13. 3[ we set the heating and coohng rate due to molec- 
ular transitions to zero in cells with ionization fractions 
between 1% and 99%. Adding the rates from local pro- 
cesses together with the photoionization rates gives the 
total rate of change of the gas energy {Q — cy-^™- and the 
neutral density (7^ - J)"'" at time i"'™. 

Once we know the rates of change of the energy and 
neutral density, we compute the time step for this itera- 
tion, 



At"-'"' = 



1 

To' 



if- PnY 



{G - ' (7^ - X)".'" ' (7^ - xy 



,(32) 



where the maximum is over the three quantities in the 
parentheses computed over all cells. Thus, the itera- 
tion time step is chosen so that the maxmium fractional 
change in the internal energy e, neutral density p„, or ion 
density p — Pn in any cell is 10%. We also want to ensure 
that the total radiation time we advance is no larger than 
the Courant time step in the initial configuration, so if 
r^" + Ar'™ > tn + Af^., we set At"'™ = t„-|-AiJ^-i"'"'. 

Once we have fixed the time step At"'™ for this iter- 
ation of the radiation update, we actually perform the 
update to find the new state 



n,m+l 









g- 



(33) 



We terminate iteration if one of the following conditions 
is satisfied: (1) we have advanced over the total time 
determined by the Courant time step in the initial con- 
figuration, i.e. t"'™+i = t" + Atg,; (2) we have ad- 
vanced over a total time that is equal to or larger than 
the Courant time step in the current configuration, i.e. 
t"'™+i > r At^l.'"; (3) the internal energy of a ceU has 
changed by more than a factor of / from the its initial 
value, i.e. max(e"'™/e", e"/e"'™) > /, where again the 
maximum is over both quantities and over all cells. If 
none of these conditions are met, we do another itera- 
tion, starting with a calculation of the photoionization 
rate by traversing the tree. 

The state when iteration terminates is the intermediate 
state we use as an input to the conservative MHD update. 
Thus, if iteration terminates after M cycles, then 



M-l 



q".* = q" + ^ At"'" 



I I ^ 



\n-i) 



(34) 



The time step for the MHD update is then set to At — 
t", the time for which we have computed the ra- 



V 



..M 



diation update. If the iteration terminates due to condi- 
tion (1), this will be the Courant time step of the initial 
state. At = Atp. If it terminates due to condition (2), 



this will be At « At^'*'^, the Courant time step of the 
intermediate, post-radiation state. If it terminates due 
to condition (3), the time step will be At = At^, the time 
it takes for the internal energy to change by a factor / in 
the cell that changes by the largest factor. In practice. 



the third condition is usually the most restrictive, since 
as the ionization front advances over a single cell the in- 
ternal energy of that cell changes by a factor of ~ 10^, 
and we find that / <C 10'^ is required to obtain good so- 
lutions. We present tests to determine what value of / 
produces acceptable results for simple radiation hydro- 
dynamics problems in § 14.3.21 

3.3. Cooling in Mixed Cells 

As part of our algorithm, we set Aki and Fki, the cool- 
ing and heating rates due to molecular processes, to zero 
in cells with ionization fractions between 1% and 99%. 
We do this in order to prevent a problem with overcool- 
ing in mixed ionization cells. To understand why this is 
necessary, let us first consider what happens at an ion- 
ization front in reality. In the fully molecular gas outside 
the front, heating due to cosmic rays and interstellar UV 
and cooling due to molecular emission balance and the 
gas is effectively isothermal; in the fully ionized gas on 
the other side, heating by ionizing photons and cooling 
by recombination and forbidden line emission also bal- 
ance to keep the gas isothermal. The ionized and molec- 
ular regions are separated by a PDR and an ionization 
front. The true thickness of the ionization front is of 
order the mean-free path of an ionizing photon, which 
in a gas at a typical molecular cloud density nn = 100 
cm^"^ is [anuy^ « 100 AU. Once an expanding ioniza- 
tion front becomes D type and sweeps up a dense shell 
of material, the thickness will be even smaller. The ion- 
ization fraction will be between 1% and 99% only within 
this thin layer. Because the hot ionized gas on one side 
of the front and the cool molecular gas on the other side 
are physically separated, and there is no mechanism to 
efficiently transport heat across the front, gas that gains 
energy from ionizing photons cannot lose it via molecular 
cooling. 

Now, however, consider what happens in a simulation 
with finite resolution. If one simulates a region that is 
parsecs or more in size, the ^; 100 AU width of an ion- 
ization front will generally be much smaller than the size 
of a computational cell. As a result, the ionization front 
will be smeared out into a width comparable to cell size, 
and this numerical mixing will supply a mechanism to 
transport heat from ionized gas into molecular gas. This 
can lead to a dramatic overestimate of the energy loss 
rate, because molecular emission is so efficient. In effect, 
the rate limiting step in cooling will become how rapidly 
energy can be transported into neutral gas by numerical 
mixing. We demonstrate this effect in a real simulation 
in§|4Xl 

Ideally one would solve this problem by either resolving 
the width of the ionization front (or at least coming close 
enough to render the overcooling region too small to do 
much damage), or using an interface-tracking technique 
to follow the position of the front on a subgrid scale. 
However, the former option is prohibitively expensive in 
thre e dimensions (though it is feasible in two dimensions, 
e.g. [Arthur fc Hoard l2006( l. and the latter is extremely 
cumbersome in a three-dimensional simulation. There- 
fore we adopt a much simpler solution: since we know 
physically that cells at an ionization front should con- 
tribute negligibly to the overall cooling rate, we simply 
set the molecular cooling rate in these cells to zero. This 
leads to a very minor underestimate of the cooling rate. 
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but the error is small because in reality the front has a 
tiny volume and therefore should contribute negligibly to 
the overall cooling rate. As we show in 14.3.41 this proce- 
dure enables us to match analytic solutions to very high 
precision, whereas if we leave the cooling on in mixed 
cells we find large errors relative to analytic solutions. 

Our method does assume that the boundary between 
the ionized and molecular gas can be described well as a 
single, simple front. If instead the molecular and ionized 
gas are mixed into a two-phase medium with a complex 
interface that has structures small compared to a compu- 
tational cell (for example dense clumps of cold molecular 
gas surrounded by a hot, ionized interclump medium), 
then our approach of simply turning off molecular cool- 
ing will fail. In this case, one would either need to resolve 
the complex interface, or use a detailed subgrid model of 
heating and cooling to account for the two-phase nature 
of the medium. 

4. ACCURACY AND CONVERGENCE TESTS 

Here we present tests of our code against analytic solu- 
tions. Since we are not aware of any analytic solutions for 
ionizing radiative transfer with magnetohydrodynamics 
in more than one dimension, all of the tests we present 
in this section are hydrodynamic rather than magneto- 
hydrodynamic, i.e. B = 0. We present magnetohydro- 
dynamic results in § [Sj 

4.1. R Type Ionization Fronts without Recombination 

A first, basic test of our implementation is the pro- 
pogation of a spherical ionization front in a medium in 
which there is no hydrodynamic evolution and no recom- 
binations. Physically, this is the type of expansion that 
should occur at the very earliest stages of an HII region's 
life, before there has been time for either hydrodynamic 
motions or for a significant number of recombinations to 
occur. In this case, a source with ionizing luminosity 
s placed in a uniform medium with density p produces 
an ionized region whose radius is simply set by equating 
the number of ionizing photons produced up to a time t 
with the number of hydrogen atoms within the radius r.^. 
Thus, we should find 



3/igs 



1/3 



t 



1/3 



(35) 



This problem is effectively a test of photon conserva- 
tion in our code. Since our radiation method is in prin- 
ciple exactly photon-conserving, we should be able to 
compute the position of the ionization front accurately 
to the resolution of our computational grid. We preform 
a test with a source of ionizing luminosity s = 4.0 x 10^^ 
ionizing photons placed in a uniform medium of den- 
sity p = 2.34 X 10^^^ g cm^^ {nn — 100 cm^^) and an 
initial temperature of 20 K. The source is at the cen- 
ter of a computational domain that runs from —10 pc 
to 10 pc in each direction. For this test we disable the 
MHD update in our code, and set the rates for all radia- 
tive processes except photoionization to zero. We use a 
resolution of 64^ cells. 

Figure [2] shows the computed radius of the ionization 
front versus time in our test. As the plot shows, the 
computed solution agrees with the analytic prediction to 
better than a percent, which is an error much smaller 




Fig. 2. — Ionization front radius versus time t (upper panel), 
and error in simulated relative to the analytic value versus time 
(lower panel) in a simulation of an R type ionization front with no 
recombinations. The upper panel shows the radius computed in 
the simulation (plus signs) and the analytically computed radius 
(solid line). We define the ionization front radius as the mean 
radius of the centers of all cells for which 0.01 < Pn/p < 0.99. The 
error is defined as Ir^i^ - r-analytl/^analyt, where r^im is the radius 
found in the simulation and Tanalyt is the analytically computed 
value. 




r (pc) 



Fig. 3. — Ratio of neutral to total density pn/p versus radius 
r for a subset of cells in our computational domain at a time 
of 10" s (3.17 kyr) in a simulation of an R type ionization 
front without recombinations. The scale bar above the location 
of the ionization front shows a ranee of ±(V3/2)Ax about the 
analytically computed radius of the ionization front at this time. 
The minimum ratio of 10"^" in the ionized region is set by a 
numerical floor in our code. 



than the size of a computational cell. Figure [3] shows 
the radial profile of the ionization fraction versus radius 
for every cell in our computational domain at a time of 
10^"'^ s (3.17 kyr), and Figure [4] shows the neutral fraction 
in an equatorial slice through the domain. As the plots 
show, the ionization front is entirely confined to cells at 
the analytically computed front radius at this time, 6.9 
pc. Thus, our code passes this test quite well. 
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Fig. 4. — Ratio of neutral to total density Pn/p in an equatorial 
slice through our 3D computational domain in a simulation of 
an _R type ionization front without recombinations. The time 
shown is the same as that in Figure [3] The dashed line shows the 
analytically-computed front position. 



4.2. R Type Ionization Fronts with Recombination 

A next test of our code is in the case where there are 
still no motions, but recombinations do occur. Physi- 
cally, this represents the next phase in the expansion of 
an R type ionization front, when the recombination rate 
becomes competitive with the photoionization rate, but 
before the excess pressure in the ionized gas causes the 
HII region to begin hydrodynamically expanding. If the 
recombination coefficient a^^^ is constant, this problem 
has an analytic solution: at time t, the radius of the 
ionization front is 



1 



-t/trec 



where 



1/3 



1/3 



is the Stromgren radius and 



tlQ 



P 



(36) 



(37) 



(38) 



is the recombination time scale for the gas. At early times 
this solution matches the case of powerlaw expansion dis- 
cussed in § 14.11 but as times t ^ tree the ionization front 
slows down due to recombinations. At times t 3> tree, the 
front radius approaches and the front velocity drops 
to zero. 

For any realistic cooling curve, for which the heaitng 
and cooling rates depend on the density, radiation flux, 
and the ionization state, there will be some small temper- 
ature variation in time and space within the HII region, 
so a^^^ will not be exactly constant and the analytic so- 
lution will not apply exactly. To enable us to compare 
with the analytic solution as well as possible, and there- 
fore to make the strongest possible evaluation of the code, 
for the purposes of this test we set a*^'^-' — 2.59 x 10~^^ 
cm~^ s^^ independent of temperature. 

We run this test using the same initial conditions and 
resolution as in § 14. 1[ which give Ts = 5.0 pc and t^cc = 
1.2 kyr for our chosen value of a^^\ For this test, we 
again disable the MHD update step in our code, but 




Fig. 5. — Ionization front radius ri versus time t (upper panel), 
and error in simulated relative to the analytic value versus 
time [lower panel) in a simulation of an R type ionization front 
with recombinations. The upper panel shows the radius computed 
in the simulation (plus signs), the analytically computed radius 
(solid line), and the analytic prediction for the radius if there are 
no recombinations (dashed line). We define the ionization front 
radius as the mean radius of the centers of all cells for which 
0.01 < p„/p < 0.99. The error is defined as - ranalytl/^'analyt, 
where rgim is the radius found in the simulation and r^nalyt is the 
analytically computed value. 



r (pc) 

Fig. 6. — Ratio of neutral to total density pn/p versus radius 
r for a subsample of cells in our computational domain at a time 
of 10-*^^ s (3.17 kyr) in a simulation of an R type ionization front 
with recombinations. The scale bar above the location of the ion- 
ization front shows a range of it(\/3/2)Aa; about the analytically 
computed radius of the ionization front at this time. 



allow the entire radiation module to operate normally. 
We run for 10^"'^ s (3.17 kyr), which is more than 2t^cc- 
Figure [5] shows the position of the ionization front versus 
time, computed as in W.ll and Figure [5] shows the radial 
profile of the neutral density at a time of 10^^ s. As we 



found in § 14.11 the ionization front is entirely confined to 
a radial extent of a single cell, so the front is as spherical 
as possible given our Cartesian grid. The accuracy of the 
front position is ~ 1%, better than a single cell. 

4.3. D Type Ionization Fronts 
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The tests described in 14. H and ? 14.21 demonstrate that 
our radiation r nodule produces correct results on its own, 
and the tests in lGardiner fc Stond (|2005ll2006[ ) show that 
the MHD module works well. However, we must also 
test how the two modules work together. We therefore 
consider the problem of the hydrodynamic evolution of 
a uniform neutral medium with initial density p„ and 
temperature r„ within which there is an ionizing source. 
This problem has a well-known solution. At early times 
there is an R type ionization front: the source ionizes the 
gas out to a radius ~ Vg, and heats it to Ti « 6400 
K, but the gas does not move and the ionization front is 
not preceeded by a shock. However, the ionized region is 
overpressured relative to the surrounding gas, and it ex- 
pands on a sound crossing time scale, tg = Vs/ci, where 
Ci « 8 X 10^ cm s^^ is the sound speed in the ionized 
gas. As it expands, it sweeps up a dense shell of neutral 
material, forming a D type ionization front. The expan- 
sion is subsonic with respect to the ionized gas, so the 
ionized region has an almost uniform density and temper- 
ature. Ionization balance requires that this density vary 



as pi (X r, 



-3/2 



so the shell is driven outward by a pres- 



sure p = piC^ ^^'^ ■ III the limit where the pressure 
in the ionized region that is driving the shell greatly ex- 
ceeds the thermal pressure in the ambient medium, and 
the neutral mass in the shell greatly exceeds the mass 
within the ionized region, conservation of momentum for 
the shell impli es that after a time t it expands to a radius 
(|Spitzeilll978h 

rs. = r.[l + -J . (39) 

The thickness of the shell relative to its radius is of order 
(cn/vsh)'^, where c„ is the isothermal sound speed in the 
neutral gas and Vsh = Ci[l + {7 /4)t/ts]~^^'^ is the velocity 
of the shell. As with the analytic solution for the R type 
ionization front, this strictly applies only if the ionized 
gas sound speed Ci is constant. For a realistic cooling 
curve, we expect this to be nearly but not perfectly true. 

We run simulations with an initial density p = 2.34 x 
10~^^ g cm~^, temperature T = 11 K, and an ionizing 
source of luminosity s = 4.0 x 10^^ s""'^. The source is 
placed at the center of a box running from —10 pc to 
10 pc in every direction. For these values of p and s, 
and assuming an ionized region temperature of roughly 
8.0 X 10^ K, we have = 0.5 pc and ts = 0.061 Myr. 
We expect that the similarity solution will apply for times 
I ^ t/ts -^i 1000. The lower limit on the time is imposed 
by the requirement that the mass in the swept-up shell 
be much larger than the mass in the interior of the HII 
region, and the upper time limit comes from the require- 
ment that the pressure in the HII region greatly exceed 
the ambient pressure, or equivalently that the thickness 
of the shell be much smaller than its radius. We run all 
our tests for 1 Myr, or 16.3^^- 

4.3.1. Comparison to the Analytic Solution 

We begin with simulations at resolutions of 64^, 128'^, 
and 256'^, using a factor / = 4 for our heating time step 
constraint (i.e. we allow the internal energy of a cell to 
change by up to a factor of 4 between hydrodynamic up- 
dates.) Figure [7] shows the radial distribution of den- 
sities, temperatures, neutral fractions, and isothermal 



sound speeds in the gas in the simulations after 1 Myr. 
As we expect, the density and temperature in the ionized 
region are almost constant. The density peaks at the po- 
sition of the overdense shell, but the shell is spread over 
a few cells, so its width decreases with resolution. This 
spreading is to be expected. The true physical thickness 
of the shell at this time should be roughly a percent of 
Vsh, smaller than the size of a single computational cell. 
However, our code has some numerical viscosity, partic- 
ularly for motions that are not aligned with the compu- 
tational grid, and this spreads shocks over ~ 3 cells. As 
Figure [8] shows, despite this spreading the front remains 
round at each resolution. 

We next compute the radius of the dense shell as a 
function of time in our simulation. We do this by identi- 
fying all the cells whose densities are at least 10% larger 
than the initial density, and then taking the average of 
their radii. Using an overdensity cutoff different from 
10% modifies the results in detail, but not qualitatively. 
We compare this measurement to the analytic solution 
in Figure [HI As the plot shows, the radius of our simu- 
lated shell matches the analytically computed radius at 
late times quite well in all the runs. The error is always 
smaller than cell in size, which is the best that can be 
expected using a Cartesian grid with finite numerical vis- 
cosity. This error is ~ 3% of the analytically computed 
radius at the lowest resolution and less than 1% at the 
highest resolution. The error generally decreases with 
time, both because the ratio of the radius to a cell size is 
increasing, and because the analytic solution is becoming 
a better and better approximation as time increases. 

4.3.2. Heatmg Time Step Convergence Tests 

We have shown that our method can match the ana- 
lytic solution for the propogation of a D type ionization 
front with a heating time step constraint of / = 4, mean- 
ing that we allow the pressure in a cell to change by at 
most a factor of 4 between hydrodynamic updates. How- 
ever, we would like to explore a range of heating time 
step constraints to determine if the accuracy of the solu- 
tion is affected by the choice of /. We therefore re-run 
our 128^ simulation with / = 2, / = 10, and / = 100 
to see if the quality of the solution changes significantly. 
The choice / = 100 corresponds to essentially no con- 
straint on the hydrodynamic time step other than the 
ordinary Courant condition. We wish to use as large a 
value of / as possible, because that will minimize the 
computational cost by enabling us to take fewer of the 
more expensive MHD updates per radiation update. Of 
course if / is too large then the cost of the radiation up- 
date will dominate the total computation cost, and there 
will be no additional advantage in increasing /. Where 
exactly this crossover occurs depends on the number of 
processors and the geometry of the ionization front. 

Figure \W\ shows the radius of the shell versus time as 
computed in our simulations with varying values of /, 
and as computed analytically. The results show that the 
radius approaches the analytic solution somewhat fall 
all /, but that the error is larger and the convergence 
slower for larger values of /. This suggests that numer- 
ical methods in which the temperature chan ge between 
hydrodynamic updates is uncon strained (e.g. lDale et al.l 
[2005tlMellema et al.l[2005l. [200l or is limit ed only to val- 
ues of / » 10 (e.g. iMac Low et al.|[2006D may produce 
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Fig. 7. — Hydrogen number density, temperature, neutral fraction, and isotliermal sound speed versus radius for a selection of 
computational cells in our simulation of a D type ionization front, at resolutions of 64^ (Ze/t column), 128^ {middle column), and 256'' 
(right column). We show only a subsample of the gridpoints to minimize clutter and file size. The points are plotted at a time of 1 Myr 
{n.7ts) after the start of the simulation. The scale bars in the uppermost panels indicate a range of it3(\/3/2)Ax about the analytically 
computed radius of the dense shell at this time, to indicate a width of 3 cells, about the expected width of a shock due to numerical viscosity. 



quantitatively incorrect results for the expansion rates of 
D type ionization fronts, although the error is likely to 
be only a few cells. This problem may not affect methods 
that implicitly upd ate the hydrodynamics and the radi- 
ation together fe.g. iMiniati fc Colellal[2006l ). rather than 
in an operator split fashion, but we are unaware of any 
such methods for ionizing radiation transport as opposed 
to local heating and cooling functions. 

Note, however, that even with large / the error does 
decrease in time. This result is easy to understand in- 
tuitively. Once the front is well into the D type phase, 
the Courant time step is set by the sound speed in the 



ionized region, which is roughly constant. As the front 
sweeps up more material and slows down, the rate at 
which gas is heated from the inital temperature to the 
ionized temperature decreases, so, on average, the frac- 
tional change in cell temperature per fixed Courant time 
step is a decreasing function of time. Thus, the error 
one makes by allowing a large change in gas temperature 
per Courant time step decreases with time. This means 
that codes with unconstrained temperature changes per 
hydrodynamic update are likely to be most inaccurate at 
early times when fronts are moving at speeds close to the 
ionized sound speed, and are most accurate at late times 
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Fig. 8. — Slice through the xy plane showing log of hydrogen 
number density in a simulation of a D type ionization front 1 Myr 
(17.7ts) after the start of the simulation, at a resolution of 64^ 
{top panel), 128^ {middle panel), and 256'' {bottom panel). The 
slices in the xz and yz planes are nearly identical. 



when fronts move slowly. 

4.3.3. Ray Rotation 

To demonstrate why it is necessary to rotate the orien- 
tation of the rays in our calculation regularly, we repeat 
our 128'^, / = 4 test calculation without ever rotating the 
rays. The mean radial positions of the dense shell and 
the ionization front are essentially the same in the two 
calculations, but the shape of the ionization front is sig- 
nificantly different. Figure fTD shows the neutral fraction 
in slices through the simulations without and with ray 
rotation after 1 Myr of evolution. With rotation of the 
rays, as shown in the lower panel, the ionization front 
boundary is smooth and round. Without rotation, as 
shown in the upper panel, there are projections of par- 




FlG. 9. — Radius of the dense shell around an HII region verus 
time {top panel), computed analytically {solid line) and by simula- 
tions at resolutions of 64^ {plus signs), 128^ {asterisks), and 256'' 
{diamonds). We also show the difference between the simulated 
and analytically-computed shell radii, r^im and r^nalyti normalized 
to the analytic radius {middle panel) and normalized to the size of 
a computational cell {bottom panel). 
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Fig. 10. — Radius of the dense shell around an HII region 
verus time {upper panel), computed analytically {solid line) and by 
simulations at a resolution of 128^, with / = 2 {plus signs), / = 4 
{asterisks), / = 10 {diamonds), and / = 100 {triangles). We also 
show the fractional error in the position of the shell relative to the 
similarity solution {lower panel), defined as |»'sim — ^analyt l/^analyti 
where rsim is the simulated radius and ranalyt is the analytically 
computed radius. 

tially neutral gas several cells into the ionized region, and 
the front boundary is much less smooth. If one is inter- 
ested in studying instabilities or similar phenomena at 
the ionization front boundary, the presence of these fin- 
gers of neutral gas is potentially problematic, since they 
may serve as seeds for instability. 
These projections are caused by differences in ioniza- 



MHD HII Region Evolution 



13 



101 — ' — " — ' — I — ' — ' — ^ — ' — ' — ' — ' — r^-^ — ' — ' — I — ^0.0 




-5 - 



-10 I , ^ , I , , . , , , , . , , ■-2.0 

-10 -5 5 10 

X (pc) 

Fig. 11. — Slice through the xy plane showing log of the ratio 
of neutral density to total density in a simulation of a D type 
ionization front 1 Myr (17.7ts) after the start of the simulation, 
with {upper panel) and withpout (lower panel) rotation of the rays. 

tion rate for cells at the same radius due to discretization 
in angle. When the orientation of the rays is constant 
over long times, small differences in ionization rate at 
different angles can build up to produce the finger-like 
structures shown in the figure. One could reduce the 
discretization error by using a more refined ray tree, for 
example by requiring that the solid angle subtended by 
a cell always be at least four times as large as the solid 
angle corresponding to a single ray, rather than twice as 
large as we require. However, this is potentially expen- 
sive, since every factor / increase in the angular reso- 
lution multiplies the time required to compute the pho- 
toionization rate by a factor of and the photoioniza- 
tion rate is computed during every radiation iteration. 
It is also overkill, since the angular disretization errors 
only build up over many time steps. Rotating the rays at 
intervals of 5 hydrodynamic time steps, as we have done 
for the simulation shown in the lower panel of Figure fTTl 
eliminates or vastly reduces the occurence of neutral gas 
fingers in the ionized region. The computational cost of 
rebuilding the ray tree every 5 hydrodynamic time steps 
is completely negligible and, for parallel calculations, re- 
quires no additional inter-processor communication. 

4.3.4. Ionization Front Cooling 

We next explore the question of cooling in cells with 
mixed ionization states. As discussed in § 13.31 we set 
the molecular cooling rate Aki = Fki = in cells where 
the ionization fraction is between 1% and 99% in order 
to prevent excessive cooling due to numerical mixing of 
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Fig. 12. — Radius of the dense shell around an HII region verus 
time (upper panel) , computed analytically {solid line) and by simu- 
lations at a resolution with cooling in mixed-ionization cell disabled 
(plus signs) and enabled (asterisks). We also show the fractional 
error in the position of the shell relative to the similarity solution 
(lower panel), defined as \rsim - J^analyt I /f analyt . where r^im is the 
simulated radius and r^^^iy^^ is the analytically computed radius. 
For comparison, we also show the fractional error corresponding 
to a difference of (v^/2)Ax (dashed line), the distance from the 
center of a computational cell to its corner. 

ionized and neutral gas. To demonstrate why this is nec- 
essary, we repeat our simulation of a D type ionization 
front with / = 4 with cooling allowed in all cells regard- 
less of ionization fraction. 

Figure [T^ shows the result. The run where all cells can 
cool lags the analytic solution even at early times, and 
the error does not decrease with time. The error at 1 
Myr of evolution is roughly 14%, which is far larger than 
a single cell; the error in the slope of the of the rgh versus t 
curve at 1 Myr is 10%. In contrast, the error with mixed- 
ionization cooling disabled is 1.5%, which corresponds 
roughly to the center-to-corner distance of a single cell. 
The slope of the curve differs from the analytic value by 
less than 1%. 

The lag is a result of numerical mixing at the ionization 
front coupled with strong molecular cooling. To demon- 
strate this, we plot in Figure [13] both the total cooling 
rate C = Axin^^ -I- Ai-cc"-e'T'i/+ + Aion-ff?^e"■^f+ (panel a) 
and the molecular cooling rate Akiti^ (panel b) as a 
function of position from our run with cooling in mixed 
cells allowed, at a time of 1 Myr. As the plot shows, with 
cooling allowed in mixed cells, the energy loss rate from 
cells at the ionization front is much larger than that from 
cells either inside or outside the front. 

One might think this is to be expected, since the cells 
at the ionization front are the ones furthest out of equilib- 
rium, and are the places where much of the energy from 
ionizing photons is deposited. However, 91% of the cool- 
ing in the region with an ionization fraction larger than 
1% comes from molecular cooling rather than from re- 
combination or forbidden line cooling, the processes that 
should dominate in ionized regions. Molecular cooling 
dominates even if we consider only much more strongly 
ionized regions. Only if we limit our attention to gas 
where the ionization fraction is greater than 74% do we 
find that molecular processes provide less than 50% of 
the total cooling rate. That molecular cooling is the 
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Fig. 13. — Slice through the xy plane showing the total cooling 
rate (panel a) and the cooling rate due to molecular processes alone 
{panel b). Note that, in order to make the cooling rate in the fully 
ionized region visible, we have selected a color scale that saturates 
in the rapidly cooling ionization front. The maximum cooling rate 
reaches > 10~^^ erg cm~^ s^^. 

dominant cooling process even in gas that is up to 74% 
ionized is a clear sign that the cooHng rate is being arti- 
ficially enhanced by numerical smearing of the ionization 
front, since for a real ionization front there should be no 
molecules to cool in regions that are significantly ionized. 

This suggests that simulations can only compute the 
correct expansion rate for HII regions if they suppress 
the effects of numerical mixing at the ionization front. 
This may explain why some numerical simulations with- 
out this precaution prod uce results that lag th e true so- 
lution significantly (e.g. iMac Low et al.l l2006l figure 1) 
unless they have enough resoluti on to at least m arginally 
resolve the ionization front (e.g. [Arthur fc Ho arc 2006), 
something that is not generally feasible in three dimen- 
sions. Note, however, that the strength of this effect 
depends on the existence of extremely rapid cooling pro- 
cesses for the non-ionized gas. If the chemistry is such 
that the neutral gas is not able to cool rapidly, the er- 
ror will be much smaller. For example, with primordial 
cosmological chemistry the lag is much smaller because 
the neutral gas cooling is much less efficient (D. Whalen, 
2006, private communication). 

An important point is that the evolution in the case 
with overcooling is qualitatively correct. There is an ex- 
panding, spherical Hll region driving a dense shell of 
swept-up material, and following a radius versus time 
curve that resembles the correct solution. It is only when 
one quantitatively checks the numerical method against 
an analytic solution that the problem becomes clear. 
This highlights the importance of performing quantita- 
tive comparisons of numerical and analytic results. 



5. HII REGION EVOLUTION IN UNIFORM MAGNETIZED 

MEDIA 

Now that we have determined the constraints our 
numerical method must obey from tests against ana- 
lytic solutions in the hydrodynamic case, in this sec- 
tion we demonstrate our method in the simplest mag- 
netohydrodynamic case. Several authors have con- 
sidred the structure of magnetized ionization fronts in 
the past. However, these treatments have been lim- 
ited to math ematically rigorous but one-dimeiisional 
analyses (e.g. iRedman et al.l ll998 HWilliams et~an i2000l: 
IWilliams k. Dvson |2001D or to sin iple, quasi-einpirica l 
models fe.g. ICarlqvist e t al.' '2003': ' kvutov et all l2005( l 
in more dimensions. [Bertoldi (1989) gives an analytic 
model for the implosion of a spherical magnetized gas 
clump exposed to an external ionizing source, but his so- 
lution is limited to that case. To date, there has been 
no detailed simulation of the evolution of an HII region 
driven by a point source expanding into a magnetized 
medium in three dimensions. 

We consider a neutral medium with initially uniform 
density pn and isothermal sound speed c„, threaded by a 
magnetic field B„ oriented in the x direction. We place 
a source of constant ionizing luminosity s at the origin 
that begins radiating at time t — Q. We restrict our- 
selves to the case most relevant to ionizing sources in 
magnetized molecular clouds, in which the thermal pres- 
sure in ionized gas at the initial density greatly exceeds 
the magnetic pressure, which in turn greatly exceeds the 
thermal pressure in neutral gas at the initial density. 
Thus, Ci 3> ^ c„, where va is the Alfven speed in 
the neutral gas and Ci is the ionized gas sound speed. 

5.1. Analytic Evaluation 

We can sketch out an approximate evolutionary sce- 
nario for this configuration using simple analytic esti- 
mates. At times f^ts, the gas has not yet had time to 
move in response to photoionization heating. Since pho- 
tons do not feel the magnetic field, the evolution should 
be identical to that of a non-magnetized R type ioniza- 
tion front. At t ~ ts, the ionized gas begins to expand 
due to its high thermal pressure. Since by assumption 
Ci ^ Va, the magnetic field in the gas is initially unable 
to resist this expansion. Thus, as in the non-magnetic 
case, a D type ionization front forms, gas moves radially 
outward, and a spherical dense shell forms. Ionization 
balance requires that when the shell bonding the HII 
region reaches a radius Tgh, the gas that remains in its 
interior and is not in the shell wall must have started 
within a distance 
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of the ionizing source. Since the magnetic field is frozen 
into the gas, the magnetic flux passing through the inte- 
rior of the ionized region is therefore 
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where i?„ is the magnitude of the magnetic field in the 
neutral gas outside the shell. In contrast, the total flux 
passing through the ionized region and the dense shell 
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combined is simply the initial flux that passed through 
a circle of radius rsh , 

= ^S„r3\ = 7ri?„r2 (^1 + ^ j ■ (42) 

Since the magnetic flux through the HII region interior 
rises as t^^'^ at late times, while the total flux through 
the HII region interior and the dense shell varies as t®/'', 
the fraction of magnetic field lines passing through the 
HII region that pass through the shell interior must de- 
cline with time as t~^/'^ . Physically, this occurs because 
magnetic field lines are being dragged with the gas out 
of the HII region interior and concentrated in the dense 
shell. 

As time progresses, the magnetic field changes the evo- 
lution in two ways. First, magnetic pressure and tension 
oppose the expansion of the shell perpendicular to the 
magnetic field, causing it to deform and become aspheri- 
cal. This effect should become significant when magnetic 
pressure is comparable to thermal pressure in the ionized 
gas, i.e. when pnv\ ^ picf. Since pi « irsh/rs)~^^'^Pn 
before magnetic effects are significant, we should see sig- 
nificant deformation of the HII region when it reaches a 
radius of order rsh ~ (ci/uA)^^'^''s, and the shell velocity 
is Vsh ~ VA- Second, magnetic pressure in the swept-up 
shell will limit compression in the shock perpendicular 
to the magnetic field to a factor of order Vgh/vA- This 
will make the shell thinner along the field and thicker 
perpendicular to it, an effect that will become noticable 
when Vsh ~ VA- Thus, we expect both magnetic effects 
to become significant when the region has expanded to a 
radius of roughly 




(43) 



where we have defined the magnetic critical radius 
for the HII region. Assuming r„i ^ r^, this happens at 
roughly a time 

^ I (^^ j U (44) 

after the start of the evolution. 

The evolution will then enter a new stage in which the 
thermal pressure inside the ionized region is small com- 
pared to the magnetic pressure, but both still greatly 
exceed the thermal pressure in the neutral gas. Since 
|Bp :s> picf, gas will be able to move only along field 
lines, and those field lines will be approximately straight. 
Since expansion is subsonic with respect to the ionized 
gas and we are concerned with times 3> tg, the density 
in the ionized gas along each field line is constant. Along 
the field, gas motions are unrestricted and the evolu- 
tion will be a normal D type ionization front as seen in 
non-magnetized gas. Perpendicular to the field, mag- 
netic pressure effects become more and more dominant 
as time goes on. The ionized gas continues to escape 
along field lines, so the pressure driving expansion per- 
pendicular to the field continues falling, which causes the 
front to slow and the density contrast to decrease perpen- 
dicular to the field. At the same time, everywhere that 
the magnetic field is at all oblique to the front, a slow 
mode will move into the neutral medium. This mode 



will remove magnetic support from the front by bend- 
ing the fiel d lines so that they b ecome perpendicular to 
the shock ()Winiams et al.ll2000D . This loss of magnetic 
support will cause the shell to transition from thick and 
magnetic pressure-dominated to thin and without much 
magnetic support over a larger and larger solid angle as 
time goes on, until only the circle exactly perpendicular 
to the initial magnetic field is not occupied by a dense 
shell. 

5.2. Simulation 

We simulate an initially neutral uniform medium with 
density p„ = 2.34 x 10"^^ g cm"^ {uh = 100 cm"^) and 
temperature T„ = 11 K, threaded by a magnetic field 
B = 14.2 ic pG (4.0 x 10~^ in the units we use in this 
paper). This gives an Alfven speed of 2.6 km s'-^ and a 
sound speed of 0.20 km s~^ in the gas at the initial time. 
The computational domain runs from —10 pc to 10 pc 
in every direction, and there is an ionizing source at the 
origin with luminosity s ~ 4.0 x 10*^ s~^, giving an initial 
Stromgren radius — 0.5 pc and sound crossing time 
ts — 0.056 Myr, assuming a sound speed Ci — 8.7 km s~-^ 
in the ionized region. For these parameters the magnetic 
critical radius is Vm — 6.0rs = 2.5 pc and the magnetic 
critical time is = QAtg = 0.53 Myr. The simulation 
uses a resolution of 256"^ cells and a heating/cooling time 
step constraint of / = 4. We run the simulation for a 
time of 3tm = 1.58 Myr. 

5.2.1. Structure at t <^ t,n 

Figures [14] and [15] show slices through the computa- 
tional domain at a time of i « tm/3 = 0.177 Myr. As 
expected, the HII region almost spherical, although there 
is already a slight asymmetry visible. The outer bound- 
ary of the dense shell is nearly perfectly round, and is 
at essentially the same radius as we would expect for 
a non-magnetized region at the same time. However, 
the shell is slightly thicker in the direction perpendicu- 
lar to the field than in the direction along the field, due 
to magnetic pressure opposing compression. Also as ex- 
pected, because magnetic field lines are being advected 
with the gas, the magnetic field inside the HII region is 
much weaker than the background field, and the mag- 
netic field in the dense shell is much stronger than the 
background field. 

5.2.2. Structure at t ^ tm 

Figures [TH] and [T7] show the same slices as Figures [Til 
and [15] at a time t tm = 0.53 Myr. As expected, 
the magnetic effects have now become pronounced. The 
inner, ionized region is strongly prolate, so that it is 
roughly a factor of 2 longer along the field than per- 
pendicular to it. Along the field, the dense shell is still 
only a few cells thick, the same as in the non-magnetic 
case, but in the directions perpendicular to the field its 
thickness is not much smaller than its radius, again as we 
expect. However, the outer edge of the expanding shell 
is still not too far from spherical, and in the direction 
along the field its radius is close to r^, the radius we 
would expect were there no field present. 

Interestingly, the shell extends somewhat further in the 
direction across the field than in the direction along it, 
so the outer radius of the swept up shell in the y and 
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Fig. 14. — Slice through the xy plane showing log of hydrogen 
number density in a simulation of an ionization front expanding 
into a magnetized medium after a time of roughly tm/3 = 0.177 
Myr. The lines are magnetic field lines. The density of lines outside 
of the HII region corresponds to a magnetic field strength of 14.2 




Fig. 15. — Slices through the yz plane showing log of hydrogen 
number density {upper panel) and log of magnetic field strength 
{lower panel) in a simulation of an ionization front expanding into 
a magnetized medium after a time of roughly tm/3 = 0.177 Myr. 

z directions is noticably larger than Tm- This is likely 
because at t = tm the expansion speed is comparable 
to the Alfven speed, and is slower than the fast mag- 
netosonic speed. Across the field, signals from the shell 
can propogate by fast magnetosonic waves more rapidly 
than the shell can expand along the field. This creates a 
density disturbance, carried by the fast mode, ahead of 
the radius that the shell would have reached had there 
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Fig. 16. — Slice through the xy plane showing log of hydrogen 
number density in a simulation of an ionization front expanding 
into a magnetized medium after a time of roughly tm = 0.53 Myr. 
The lines are magnetic field lines. The density of lines outside of 
the HII region corresponds to a magnetic field strength of 14.2 fiG. 




Fig. 17. — Slices through the yz plane showing log of hydrogen 
number density {upper panel) and log of magnetic field strength 
{lower panel) in a simulation of an ionization front expanding into 
a magnetized medium after a time of roughly tm = 0.53 Myr. 

been no magnetic field. Since fast waves cannot pro- 
pogate parallel to the magnetic field, the leading edge of 
the disturbance has advanced less in the x direction than 
in the y or z directions. 

5.2.3. Structure at t^ tm 

Figures [T51 and [TOl show the simulation at our final time 
slice, t = 1.58 Myr = 3t„i. In the evolution after time 
tm, the inner ionized region becomes more prolate, ex- 
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Fig. 18. — Slice through the xy plane showing log of hydrogen 
number density in a simulation of an ionization front expanding 
into a magnetized medium after 1.58 Myr (3tm). The lines are 
magnetic field lines. The density of lines outside of the HII region 
corresponds to a magnetic field strength of 14.2 ^G. 



panding more rapidly along the field than perpendicular 
to it. Along the field, in the x direction, there continues 
to be dense shell of swept up material. 

In the y and z directions, perpendicular to the field 
lines, the thick shell of gas bounded between the ioniza- 
tion front and the fast mode front, has also continued 
to expand. However, its overdensity is decreasing, and 
at this time is only 120 cm~^, so it is within 20% of 
the background density, and smaller than the density at 
time tm- The shell expansion velocity perpendicular to 
the field lines is much smaller than the Alfven speed in 
the neutral gas, and is comparable to the sound speed, so 
the leading edge of the shell is no longer a shock. Instead, 
the shell is just a small overdensity that is gradually re- 
laxing away. 

The boundary between the neutral and ionized regions 
in the y and z directions appears to be somewhat cor- 
rugated, but this may well be a numerical artifact, par- 
ticularly since Figure [19] shows that the structure at the 
ionization front is aligned with the computational grid. It 
may be that for ionization fronts that are nearly static, 
as this one is due to magnetic confinement, rebuilding 
the ray tree once every five time steps is insufficient to 
prevent the growth of structures. 

In between the x direction and the y direction, there is 
a region where the field is o blique relative t o the i oniza- 
tion front. As predicted bv IWilliams et al.] (|2000l ). this 
region shows two types of behavior depending on how 
close the field is to being parallel or perpendicular to the 
front. Closer to the x direction, the field been bent so it 
is close to perpendicular to the front, and a dense shell 
bounded has formed at the ionization front. Closer to 
the y direction, the field is not bent and enters the front 
at an oblique angle, and there is no shock or dense shell. 
The solid angle over which a dense shell exists appears 
to be increasing in time. 

While expansion of the ionized region perpendicular 
to the field has slowed, the expansion velocity along the 
field lines is actually larger or about the same as it was at 
time tra- To understand the origin of this effect, consider 
the jump conditions that govern the expansion of the 
dense shell in the x direction at late times, when the 
expansion is strongly subsonic with respect to the ionized 




Fig. 19. — Slices through the yz plane showing log of hydrogen 
number density {upper panel) and log of magnetic field strength 
{lower panel) in a simulation of an ionization front expanding into 
a magnetized medium after 1.58 Myr (3tm). 

gas, but the thermal pressure in the ionized gas is still 
much greater than the thermal pressure in the neutral 
gas. In the rest frame of the shell, this means that ram 
pressure of neutral gas flowing into the shell must balance 
the thermal pressure of ionized gas downstream from it. 
Thus, PnVsh ^ Picf , SO the shell expansion velocity obeys 
Vsh fx Ciy^pi/ pn ■ Thus far our analysis could apply just 
as well to a non-magnetized HII region. The difference 
that the magnetic field makes is that it greatly reduces 
the expansion rate of the ionized region in two directions. 
This prevents the density from dropping as fast as in the 
hydrodynamic case, which keeps the expansion velocity 
larger at later times. 

However, while the decline in ionized gas density and 
pressure is slower than in the hydrodynamic case, since 
gas continues to flow out of the ionzied region along field 
lines into the end caps, the pressure inside the HII region 
is always decreasing after time tm ■ By 3 tm , the HII 
region is in the limit where the magnetic pressure inside 
the ionization front is higher than the thermal pressure. 
As a result, the field lines have begun to straighten out 
inside the ionized region, so they are kinked only at the 
dense end caps. The magnetic energy density inside the 
ionized region is actually larger than at time tm, since for 
t ?J t„i field lines that were pushed into the dense shell 
start to move back into the ionized region interior. 

6. DISCUSSION AND CONCLUSION 

We have demonstrated an algorithm for computing the 
evolution of magnetized molecular gases subjected to in- 



18 



Krumholz, Stone, & Gardiner 



tcrnal sources of ionizing radiation, which is potentially 
applicable to molecular clouds. In testing our algorithm, 
we have discovered three conditions that are likely to ap- 
ply to ionizing radiation hydrodynamic and magnetohy- 
drodynamic codes in general. First, to achieve maximum 
accuracy the update time step must be limited so that 
the temperature in cells docs not change by more than a 
factor of / < 10 between hydrodynamic or magnetohy- 
drodynamic updates. Larger values of / produce small 
but significant errors in the expansion rates of D type 
ionization fronts. Second, when using a ray-tracing ap- 
proach to compute the ionizing radiative transfer, one 
should rotate the orientation of the rays periodically to 
avoid a build-up of errors caused by the discretization of 
angles around the ionizing source. Failure to obey this 
conditions results in fronts that should be spherical devel- 
oping asphcrical features, and in more complex calcula- 
tions this could potentially seed instabilities. Third, and 
most significantly, one must avoid overcooling caused by 
numerical smearing of the ionization front. This can be 
handled most easily by suppressing cooling in cells with 
mixed ionization fractions. Failure to obey this condi- 
tion leads to an unphysical loss of energy from expanding 
HII regions that causes them to lag analytic solutions by 
tens of percent. A calculation that satisfies these three 
constraints can reproduce the analytic solution for the 
expansion of a D type ionization front to an accuracy of 
a percent for at least ^ 20 ionized sound-crossing times. 

Using our algorithm, wc report the first three- 
dimensional simulations of the expansion of an HII re- 
gion into a magnetized gas. We show that the presence 
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